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MEMORANDUM 12-4-58A 


THREE-DIMENSIONAL ORBITS OF EARTH SATELLITES, 

INCLUDING EFFECTS OF EARTH OBLATENESS 

AND ATMOSPHERIC ROTATION 

By Jack N. Nielsen, Frederick K. Goodwin, 
and William A. Mersman 

SUMMARY 


The principal purpose of the present paper is to present sets of 
equations which may be used for calculating complete trajectories of 
earth satellites from outer space to the ground under the influence of 
air drag and gravity, including oblateness effects, and to apply these 
to several examples of entry trajectories starting from a circular orbit* 

Equations of motion, based on an "instantaneous ellipse" technique, 
with polar angle as independent variable, were found suitable for auto- 
matic computation of orbits in which the trajectory consists of a number 
of revolutions. This method is suitable as long as the trajectory does 
not become nearly vertical. In the terminal phase of the trajectories, 
which are nearly vertical, equations of motion in spherical polar coordi- 
nates with time as the independent variable were found to be more suitable. 

In the first illustrative example the effects of the oblateness 
component of the earth’s gravitational field and of atmospheric rotation 
were studied for equatorial orbits. The satellites were launched into 
circular orbits at a height of 120 miles, an altitude sufficiently high 
that a number of revolutions could be studied. The importance of the 
oblateness component of the earth's gravitational field is shown by the 
fact that a satellite launched at circular orbital speed, neglecting 
oblateness, has a perigee some 67,000 feet lower when oblateness forces 
are included in the equations of motion than when they are not included. 
Also, the loss in altitude per revolution is double that of a satellite 
following an orbit not subject to oblateness. The effect of atmospheric 
rotation on the loss of altitude per revolution was small. As might be 
surmised, the regression of the line of nodes as predicted by celestial 
mechanics is unchanged when drag is included. It is clear that the 
inclination of the orbital plane to the equator will be relatively 
unaffected by drag for no atmospheric rotation since the drag lies in 
the orbital plane in this case. With the inclusion of atmospheric 
rotation it was found that the inclination of the plane changed about 
one-millionth of a radian per revolution. Thus the prediction of the 
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position of the orbital plane of an earth satellite is not complicated 
by the introduction of drag. The line of apsides, which without drag 
but with oblateness moves slowly in space, tenc.s to move with the satel- 
lite when drag is included in the calculations,: As a results, the usual 

linearized solutions based on oblateness alone must be basically altered 
when drag is included to take into account the rapid movement of the line 
of apsides. 

In the second illustrative example the final revolution was calculated 
to impact for a number of trajectories in an orbital plane inclined at 65° 
to the equator. Of particular interest is the large effect the oblateness 
gravitational field and atmospheric rotation can have on the impact point. 
For a value of CpA/m of unity, and for an initial downward angle at 
80 miles altitude of 0.01 radian, such as might be utilized for manned 
re-entry, oblateness had an influence of about 300 miles in the impact 
point, and atmospheric rotation had about a 150-mile influence. 

It was found that two-dimensional solutions neglecting atmospheric 
rotation can be used to approximate three-dimensional solutions including 
atmospheric rotation. In this connection two-c.imensional theories must 
be Interpreted as being viewed by an observer on a rotating earth. 

The importance of various terms in the equations of radial and 
tangential motion is examined for various calculated trajectories. The 
validity of the principal assumption in the approximate equation of motion 
of TN 4-276 was thus confirmed for a satellite speed less than about 
99“P e rcent circular satellite velocity. Certain gaps in our theoretical 
and experimental knowledge are pointed out insofar as they influence our 
ability to calculate complete trajectories from launch to impact. 


INTRODUCTION 


Much interest exists in the dynamics of earth satellites, and a 
number of papers (e.g., refs. 1 to 6) in the field have recently appeared. 
These papers consider only parts of the total trajectory; they are also 
usually limited to two-dimensional trajectories,, or they neglect air 
drag. The principal purpose of this paper Is to present sets of equations 
which may be used for the calculations of complete trajectories from outer 
space to the ground under the influences of air drag and gravity, including 
oblateness effects, and to apply these to several Illustrative examples of 
entry trajectories starting from a circular orbit. This purpose cannot be 
achieved at present on the basis of analytical solutions except possibly 
by patching together such solutions. Therefore, automatic computing 
machinery was used in the study. Forms of the equations of motion suit- 
able for automatic computation are presented, and a number of calculative 
examples are carried out. The calculative examples are chosen to illustrate 
the relative importance of various physical forces acting on earth 
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satellites. Some of the essential physical features of the calculated 
solutions are discussed in the hope that the resulting insight will lead 
to more inclusive analytical solutions. In connection with the numerical 
calculations the authors wish to acknowledge the contributions of Miss 
Marcelline Chartz in programing the equations for the digital computing 
machines and in supervising the computations. 


SYMBOLS 


a 

& r ,a 'K ,a '^ 

A 

— > 

A 

b 

CD 

D 

Dq),D a J 


equatorial radius of earth 

-> — ► -► 

components of acceleration vector along i r ,i^, 3X1(1 ^ 
directions 

reference area of satellite 
acceleration vector of satellite 
polar radius of earth 
drag coefficient 
drag force 

components of drag along i r ,i-^ ,i^ ,1^, and i a directions 


e 

E 

E 


eccentricity of earth 

eccentricity of instantaneous ellipse, 
total energy of satellite 


( l 2 - k 2 ' 

'v z 2 . 


1/2 


F r> F V F \|/' 


g 


G r> G V G T|r 

G cp> G a 


vector force acting on satellite 


components of F along i r and i a directions 
component of gravity acceleration at earth* s surface^ no 


oblateness 

vector force per unit mass acting on satellite because of 
gravity of earth 


components of G along i r ji^ji^jiq). 


and i a directions 


Oblate ness 

3 

Yes 

Yes 

No 

uu otuue 
atmosphere 
within range 
of oblate - 
ness forces 

2 

Circulari- 





During circu- 

T •? r. r,+- 

•a 
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vector angular velocity of i^i^,!^ system in inertial 
framework 

angular velocity of earth about polar axis 

A ft change in angular velocity due to motion of orbital plane 

— ► 

Aftf angular velocity for fixed orbital plane 
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The angle |3 is the bearing angle of the orbital plane in the r, A, 
system. The bearing angle for an earth observer is different from P 
because of earth rotation. 

The relationships among the various angles can be readily established. 
The following six relationships are useful 


sin ^ = 
sin (A - 0 ) = 

cos (A - 0 ) = 

sin a cos cp = 
cos a = 

tan a = 

The direction cosines between t 
system are 


sin Cp sm a 

sin cp cos a 
cos \|r 

cos cp 
cos \jr 

cos \|r cos p 
cos sin p 
tan 

sin(A - 0 ) 

Le r, cp, a system and the r. A, \|r 



r 
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a 

r 

1 

0 

0 

A 

0 

sin p 

-cos P 

♦ 

0 

cos p 

sin p 


Properties of the Atmosphere and the Earth 


The forces due to air drag depend on the density and rotation of 
the atmosphere, and the forces due to gravitation depend on the size, 
shape, and density distribution of the earth. The atmospheric properties 
that are adopted for the exact calculations of this report are those con- 
tained in reference 9* However, as a result cf observation of the satel- 
lite 1957 cb, Snutnik I, the ARDC densities atove about 200 KM appear to 
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will be a periodic variation in density for a purely circular orbit, due 
to the nonspherical figure of the earth, except for equatorial orbits. 

For a polar orbit, for instance, there would be density changes between 
the equator and the poles corresponding to a 13-mile change in altitude 
for a circular orbit. We consider this density change to be due to the 
nonspherical figure of the earth and do not call it an oblateness effect, 
a term applying only to the gravitational field. 


In the calculations it is necessary to know the size and shape of 
the earth. It is assumed that the earth is an ellipsoid (oblate spheroid) 
with the dimensions adopted at the International Geodetic and Geophysical 
Union of 1924 , as given in reference 8. The equator is assumed to be a 
perfect circle of radius a, and the polar radius is b. 


a = 6,378,388 ±18 meters 
= 3963*3386 statute miles 
= 20 , 926,428 feet 
b = 6 , 356 , 911 * 9^-6 meters 
= 3949 . 9941 statute miles 

= 20 , 855,969 feet 


a - b 1 
a - 297 


0 . 003367003 ^ 


( 5 ) 


e 


2 



0.0067226700 


(1 meter = 3.28083333 feet) 


( 6 ) 


It is also necessary to specify a value of the earth's gravitational 
field strength for the purposes of the calculations. If F is the gravi- 
tation attraction between two masses M and m a distance r apart, then 
the universal gravitational constant K is given by 



If g is the force on a unit mass due to the earth's gravitational field 
at the surface of the earth, M the mass of the earth, and r is equal 
to R, the "mean" radius of the earth, then 

KM = gR 2 

It turns out that we shall have use of the quantity KM which is given 
in reference 11 as 
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KM = (l ±8xl0" 5 )(3.986329xl0 2C cm 3 /sec 2 ) 
KM = 0.1 i K)77500xl0 17 ft 3 /sec 2 


(7) 


Another quantity of the earth’s gravitational field with which we 
will be concerned is its quadripole moment. Because of the nonspherical 



Sketch (a) 


mass distribution of the earth, its gravita- 
tional field is not spherically symmetric, 
but has a quadripole moment. The potential 
is that due to a quadripole formed as shown 
in the sketch. A quartet of sources and 
sinks with strengths of equal magnitude are 
placed on the polar axis of the earth near 
its center. The sources are located at the 
center with the sinks equally spaced above 
and below the center. The sinks are permit- 
ted to approach the sources while the pro- 
duct of strength and spacing remains constant. 
The net result is a quadripole of the present 
sort which resembles a dipole pair with mirror 
symmetry. Also shown in the sketch are the 
radial and tangential components of the 
gravitational field due to oblateness. The 
gravitational potential including the non- 
spherical component is usually taken as 


0 




cos 



(8) 


The force per unit mass in any direction is the gradient of 0. The 
value is the geocentric latitude and \i is a dimensionless constant 

given in reference 2 as 

6|_t = (1.637 ±o.oo4)xic -3 

From reference 11, a value of 6[x calculated for the international 
ellipsoid is 

6\i = 1.638x10" 3 

We will use this value in the calculations. 

It is probably interesting to note that the surfaces for which <£> 
is a constant are not oblate spheroids nor does one coincide with the 
assumed geometric figure of the earth. The surfaces, $ equal to a 
constant, in principle are everywhere normal to the earth's gravitational 
field. If the vector gravitational field is corrected for the accelera- 
tion due to rotation of the earth, then the resultant vector field should 
be normal to the ’’mean" surface of the oceans to prevent their flowing 
toward the equator. 


1 1. c Qa 
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If we designate^the gravitational force per unit mass in the 
i r , iq,, and i a directions hy G r /m, G^/m, G^/m, G^/m, and 

we have in the r. A, \|r system 


Gr _ ^0 _ KM 3 pKMa 2 

m Sr r 2 r 4 

Ga _ 1 S£ _ 

m r cos \|r SA 


(1-3 cos 2\|r) 




> 


G\|/ 1 So 

m ~ r St 


6KM|ia 2 

r* 


sin 2\|r 


J 


( 9 ) 


In the r, cp, a system we have 

G r KM 6pKMa 2 


m 


(l~3 sin 2 cp sin 2 a) 


Gm 12 KM|ia 2 „ 

— = - sin cp cos cp sin a 

m r 4 y r 


Ga 

m 


~4 

12 KMpa 2 


( 10 ) 


sin cp cos a sin a 


The equations are obtained, by noting that Gcp and G a can easily be 
obtained from and G^ by a rotation involving P (see fig. 1(a)). 

The quantities \|r and p are then eliminated in favor of cp and a 
through relationships derivable on the basis of spherical trigonometry. 


Air Drag Forces 


Besides the forces on the satellite due to the gravitational field 
of the earth, we will also be concerned with the air drag on the satel- 
lite with the atmosphere rotating with the earth. In the r, A, \|r 
system the velocity of the satellite relative to the rotating atmosphere, 


V R = i r r + i^r cos \|r (A - ) + i^rijr 


(in 


where Q e is the rotational speed of the earth. Since the drag is in 
opposition to the motion 


D 



m 


m 


(12) 



Ik 


The components of the drag in the r, A, t system are thus 



In the r, Cp , a system the velocity V R is (as we will show) 

V R = i r r + ifpVq, - i^ra e cos ^ 

— ► — ► — ► — ► 

= ipr+i^Vq, - ra e cos ^(i^sin p - i a cos p) 

= i r r + - rn e cos a) + i a rft (; sin a cos cp (l4) 


where we have made use of the relationships 

sin 3 cos \jr = cos a 
cos 3 cos ^ = sin a cos T 
The drag components per unit mass in the r, cp, a system are thus 

D r 1 . ASjjff 

- = - 2 pVRr \~sr 


( 15 ) 


= - \ pVr(Vcp - rft e cos a)( 


D, 


a 

m 


1 /C pA\ 

2 pV R rfi e sxn a cos cp ( — J 


(16) 


SA TELL ITE KINEMATICS AND EQUATIONS OF MOTION 
Equations of Motion in Geographical Coordinates 


The geographical coordinates r, A, ^ are a special set of spherical 
polar coordinates useful for certain problems of trajectory calculations. 
They are related to the inertial coordinates X, Y, Z as follows: 



5L ; 
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X = r cos cos A 
Y = r cos \|r sin A 


Z = r sin \)/ 


The linear velocities in the r, A, \|f directions (fig. 1(a)) are 

V r = r 

= rA cos \|r 

= r i 

and the corresponding components of the angular velocity of the 
i r , i^, i^ coordinate system are 

fir = A sin \j/ 

n A = -t > (17) 

fi^ = A cos \|/ I 

The angular velocities referred to here concern rotations of the 

^r> ^A> system in the inertial framework. The acceleration components 

are 

•• : 2 *2 2, I 

a r = r - ry -rA cos \|/ I 


A 


r cos dt 


7PF (^ 2 cos2^a) 


(18) 


a^ = r\|r + 2rijf + rA sin \|r cos T|r 


The forces per unit mass in r. A, and \]/ directions are equal to 
the accelerations in those directions, so that the equations of motion 
are : 


r r •• • 2 • 2 2 

— = r - r\[r - r A cos \[r 

^A _ 1 1 / 2 2 1 \ 

m ” r cos dt ' r cos 

F \|7 •• • . *2 

— = r\|r + 2r\|f + rA sin \|r cos \|r 

Under the actions of gravity, including oblateness, 
forces per unit mass are 


) ( 19 ) 


j 

and of air drag the 




Equations of Motion in Orbital Plane; Motion of Plane 


For certain kinds of calculations the equations^of motion ir^ terms 
of r, 9, and a are convenient. The unit vectors i r , iq), and i a form 
a right-handed system in that order. Because the r, cp, a coordinate 
system is not an inertial system, account must be taken of its rotation 
in establishing the equations of motion. It is clear that r, cp, a, 
and 0 are necessary to describe the position of the satellite, so 
that 0 will also enter the equations of motion. Let us first establish 
the kinematic relationships of the. system. Since there are four coordi- 
nates specifying the satellite position, rather than the usual three, we 
can anticipate an extra kinematic relationship involving r, cp, a, and 0. 


A convenient method for deriving certain kinematic relationships is 
to obtain the quantity ftxr by two methods^and then to equate its com- 
ponents. One convenient way to establish ftxr^ is to consider the follow- 
ing vector transformation between any vector, r, and its time rate of 
change 




( 21 ) 


as discussed, for instance in reference 12. In this equation dr/dt is 
the total rate of change of the vector, r including both changes in the 
magnitudes and directions of its components. The quantity Sr/St refers 
to the rate of change of r due solely to the changes in the magnitudes 
of its components but with no changes in their directions. The angular 
velocity, ft, is that of the axis system in which the components are 
expressed. (The ability to hold the direction of the components fixed 
presupposes the knowledge of an inertial systen to which ft can be 
referred. ) 


Let us now take r to be the radius vector of the satellite so 
that dr/dt is the satellite velocity V From the definition of the 
orbital plane , we have in the r, cp, a system 



12-4 - 5 &A 
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The quantity 8 r/ 8 t is the satellite velocity if the i r , iqj , and i a 
directions are fixed. 


8 r _ 

5t = 


i r r 


From equations (21), (22), and ( 23 ), 


we thus obtain 


(23) 


Qxr - i(pVq) (24) 

— ► — > — ► 

The second method used to obtain Qxr is to construct ft and then 
to take its cross product with r. To establish the angular velocity we 
note that of the independent variables , r, cp, a, 0 , changes in all 
quantities except r cause angular velocity of the moving axes. If we 
hold the orbital plane fixed in space , we get an angular velocity vector 
along i a due to cp. 



(25) 


Now, because of motion of the orbital plane specified by Q and a, we 
have the additional angular velocity Aft which is clearly 

Aft = ± Z Q + i x t<x (26) 

The vector Aft can easily be transferred to the r , cp, a system by 
means of the following table of direction cosines 


X* 

Y' 

Z 

Aft = i r (9 sincp sin a + a cos cp) + ^(0 cos cp sin a-a sin cp) + i a (0 cos a) 

(27) 

The total angular velocity is now 

ft = A ftp + Aft = i r (0 sin cp sin a + a cos cp) + 

iq )(0 cos cp sin a-a sin cp) + i a (<p +0 cos a) ( 28 ) 
— ► — ► 

If we take the vector product ftxr from equation (28), we get 
ftxr = -i a r (0 cos cp sin a-a sin cp) + icpi'(q) + Q cos a) 


r cp a 


cos cp 

-sin cp 

0 

sin cp cos a 

cos cp cos a 

-sin a 

sin cp sin a 

cos cp sin a 

cos a 


( 29 ) 
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Comparison of equations (24) and (29) yields a pair of relationships 


Q cos cp sin a - a sin cp = 0 


• ♦ 


cp + 9 cos a = 


The angular velocity of the satellite is thus 


Vcp 


(30) 


-► -► . sm a -> / v ct: 

a = ±yB — — - + V — 
~ sm cp 


* a - /V cp N 

t Ini ' 

r cos cp °\r j 


(31) 


The acceleration in r, cp, a coordinates is derived from a vector 
transformation equation similar to equation (21 ) 


dV _ 8V 


+ SKV 


dt eft 

— ► . — *■ ~ > — ► 

The quantity, dV/dt, with i r , i^, and i a fi>ed is 

c>V - i r + i V 

5t " X r r + V<P 

since the rotation of the moving coordinates do not contribute to dv/dt. 
The vector product QxV is 


nxv = 


icp 

0 

Vcp 


“■a 


V^ 


0 


V T rVcp 

- i r r + 1 cp r + a>‘r v cp 


Tlius the acceleration i£ 


- dV VqA -+ /. . Vm- 

A = = i]4 : r - — — J + i<pM Vy + r — J+ laVcp^r 


(32) 


The equation of motion then can he written 
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Fr _ 

r -ZC 

m 

r 

ii 

(v^) 

II 

Vq)f2 r 


subject to the kinematic equations 


• a % 

cp + G cos a = — j- 
r 


a r = 


a 


cos cp 


9 sin a 
sin cp 


(33) 


(34) 


The components of the forces due to gravity including oblateness and 
drag in equation (33) are 


F r KM 6pKM: 


m 


4 (1-3 sin 2 cp sin 2 <x) - ~ pV R r 


Fcp 12KMp.a 2 _ . 0 1 YCnA\ 

— = Bin cp cos cp sm 2 a - 5 P V R (V cp - riJ e cos a )(- -%-J 


>(35) 


F a 12KMpa 2 

~ = - a sin cp cos a sin a 

m 1 

It is noted that equations (33) and (34) are the equivalent of a set of 
six first-order differential equations which fully determine the history 
of six orbital elements for given initial conditions. 


1 / YcpA 

2 pV R (rH e sin a cos cp)' 


Equations of Motion in Terms of a Particular 
Set of Orbital Elements 


There are many sets of six orbital elements, the history of which 
in space and time specify the path of the satellite. The set to be used 
depends on the problem to be studied. For instance, during circulariza- 
tion of the elliptical orbit , it would be reasonable to use among the 
orbital elements those physical quantities such as eccentricity or length 
of the major axis which describes the ellipse. In fact, the concept of 
the instantaneous ellipse to fit the trajectory at every point has been 
used for a long time in celestial mechanics. We will now consider the 
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equation > of motion for a special 
set of orbital elements based on 
an instantaneous-ellipse technique 
This set of equations is used for 
certain trajectories presented 
herein. 

The instantaneous ellipse 
lying in the plane of the orbit is 
shown in sketch (b). The ellipse 
is characterized by the angle 5, 
the eccentricity E, and the semi- 
major axis l . The equation of 
the instantaneous ellipse in polar 
coordina:es is 


1 1 + E cos (cp - b) 

r ~ l (1 - E 2 ) 


One of the parameters of interest is the rate at which the radius 
vector of the elliptical path is sweeping out area. Since the rate is 
constant for a central force field, it should be a slowly changing func- 
tion of time when drag or oblateness is included. The rate at which area 
is swept out is £/ 2, where £ is given by 


I = rV 9 


Let us now introduce two new variables specifying the value of £ and E 
for the instantaneous ellipse which is to have the same value of £ as 
the trajectory 


P = 



q = pE 


( 37 ) 


Here | Q 
constant 


is some reference value of | . 


Lo 



Let us also introduce another 


( 38 ) 


We now wish to relate p and q to the constants of the instantaneous 
ellipse so that we can express r in terms of p and q. In "fitting" the 
instantaneous ellipse we have taken the trajectory and tlje ellipse at 
point P to have^the same values of the^radius vector, r, and of the 
velocity vector, V. Since the vectors r and V totally determine the 
dynamical state of a particle, the dynamical states are matched precisely. 
It follows, therefore, that the momentum and energy are also matched. 

The rate of sweeping for the ellipse %j 2 is constant and its value is 
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fixed by conditions on the trajectory at P. Also, the same can be said 
for the sum of the potential and kinetic energies . Now the rate of 
sweeping out of area by the ellipse is 

1 _ 

2 “ T 

where k is the semiminor axis and T is the period. The period depends 
only on the length of the major axis 



T = -p= Z 3/2 (4o) 

n/KM 

The relationship of p to the parameters of the instantaneous ellipse 
from equations (37), (38), (39), and (4o) is 


_P_ = 1 

L 0 Z (l - E 2 ) 

It is clear that equation (36) can now be expressed 

1 _ p + q cos T) 
r L 0 


(41) 


(42) 


These preceding remarks are adequate for establishing the set of 
variables that will now be used in the six equations of first order to 
describe the motion of the satellite. First we must decide on the 
independent variable. The variable cp is convenient if we are concerned 
with how quantities vary per revolution. Time as the independent variable 
is convenient for many other problems. With cp as the independent 
variable ^ the dependent variables are taken to be 


Parameters of 
instantaneous ellipse 


p = ^o 2 A 2 

q = pE 
5 = cp - 


a 

0 


Parameters specifying 
orientation of orbital 
plane 


t 
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Equations (33) and ( 34 ) consist of a system equivalent to six first-order 
equations . In terms of the six preceding variables, the equations of 
motion given here without proof are 




dt = 

1 


dcp 

U) 


dp _ 

_ _2p 

( Y A 

dcp 

|wu 

U j 

dq = 


SOS T) + S 

dcp 

dcp 


d6 _ 

d e 


dcp 

_q - dcp 

cos a - 

da 

/Fa \ ; 

cos cp 

dcp 

\ m / 

wVq) 

d G 

(l a\ . 

sin cp 

dcp 

1 m J 1 

wV^sin a 


&P . 
tt si: 
dcp 




( 43 ) 


J 


where the quantities co, Vq>, u, S, and tj folio/ from the variables and 
from the reference quantities £ 0 and L 0 : 


gj = |u 2 - ( — 


F a \ sin cp cos a 


m 


£u 2 


|u sin a 1+ (d0/dcp)cos a 


>2 b o 
; " P 

1 p + q cos rj 

U = - = ; 


1 dp 

S = q sin *n + 

2 p dcp 


n = cp - 6 


v cp r 


| 2 u 2 \m 


+ P 


x d0 

1 + — cos a 

\ dcp 


> ( 44 ) 


J 


Certain other identities arising in the developnent are also useful, 
especially in specifying the initial conditions: 
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E = - 


£ 2 p(l-E 2 ) 


. £q sin ri 
r = L 

1 2 

V 2 = r — ^ (p 2 + 2pq cos T) + q 2 ) 


p(l-E 2 ) J 

During the course of the present work there appeared an English 
translation of a Russian paper by Taratynova, reference 1, containing 
equations of motion equivalent to the foregoing set. Taratynova r s equa- 
tions were taken from an astronomy text in Russian of 1936 vintage. His 
equations differ from the present ones in several particulars. First , 
the latus rectum and the eccentricity of the instantaneous ellipse are 
used instead of p and q as dependent variables. As independent variables 
both t and t] are used. The independent variable t can be changed to 
cp, as in the above set of equations , simply by multiplication by co. 

It is Important to inspect equation (43) for singularities arising 
from zeros in the denominators. Inspection shows that such zeros can 
conceivably arise through w, Vq), a, or q. From equation (44) let us 


"[l + (d0/dcp)cos a] 


It Is clear that oo and Vq) both approach zero together if d0/dcp is 
small , the only case that need concern us. If the path of the satellite 
is vertical in the orbital plane, go and are zero. Equation (43) is 
then ill-conditioned, in the large changes in the dependent variables 
are accompanied by small or zero changes in the independent variable. 

In this instance it becomes necessary to change to a new Independent 
variable such as the time. When a approaches zero, the orbital plane 
is approaching an equatorial plane. However, the zero in the denominator 
of the equation for d^/dcp is only apparent. This is the case because 
F a in the numerator, which is due to drag or the gravitational field, 
is proportional to sin a as shown by equations (10) and (l6). 

The zero in the denominator of db/dcp due to q is of importance 
for circular orbits, as can be seen by the following expression for the 
eccentricity 

E = i = ^!i! 

P to 2 


(^ 7 ) 
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It is clear that q will be zero for circular orbits (E = 0) since V<p 
and r are not zero. To avoid this difficulty, consider two new dependent 
variables to replace q 


v = q cos 5 
w = q sin 5 


(48) 


Then we have 


q cos tj = q cos(cp -5) = v cos cp + w sin cp 


(49) 


The differential equation for q is then replaced by the following set 


with 


dv d0 dp — 

— = v — cos a - — cos cp + S sm cp 
dcp dcp dcp ^ 1 


dw 

dcp 


d 9 dp 

-v — cos a - — sm cp - d cos cp 

dcp dcp Y Y 


(50) 


_1_ |^pY 

2p \dcpy v 


sm cp 


W COS cp) + —0^2 (~ + P T"*" + COS <X ) (5l) 

U 2 u 2 V m y J V dc P / 


It is noted that the zero in the denominator of d&/dcp has now been 
eliminated. 


The equations of Taratynova also contain a singularity. How the 
singularity was handled in obtaining the tabulated results given in the 
paper for circular orbits is not discussed. 


It is clear that the system of six equalities included in equa- 
tion (43) can readily be converted from cp as the independent variable 
to t as the independent variable. It is sufficient to invert the first 
equality and to multiply the next five equalities by u>. 

It is possible to write the equations of notion in a form which is 
linear in the applied forces. In this case we also introduce the 
eccentricity, E, rather than q as dependent variable. 
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da cos cp 
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^cp 


a 
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1 

Ep 
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Ep 

2 
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~E sin 2 T] 
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pr^ sin tj 

I 
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dcp Vcp sin cp cos a ^F a 
dt ~ r £u sin a 


m 


£S> 

m 


t 2 

sp 

KM 


(52) 




The quantity in the square brackets represents the applied radial force 
less the spherical component of the gravitational field. The foregoing 
form of the equations of motion is convenient for deducing secular 
trends . 


ILLUSTRATIVE EXAMPLES 


The equations of motion as developed in terms of orbital elements 
and in geographical coordinates have their own particular uses. Two 
such uses will be illustrated in the present paper. In the first example 
we will examine the effect of the earth’s oblateness and of the air drag, 
including atmospheric rotation, on the orbits of satellites in the equa- 
torial plane. Such an example will use the equations of motion in terms 
of orbital elements. The second example in terms of geographic coordi- 
nates is concerned with the effect on the impact point of the earth’s 
oblateness and of atmospheric rotation for nonequatorial orbits. 


First Illustrative Example 


As a first illustrative example, consider a satellite of given 
CDA/m launched horizontally eastward on an equatorial orbit. Let the 
initial height about the earth be given, and let the velocity be that 
for a circular orbit without consideration of oblateness effects 



2 6 


(Vo = 


(V C ) 0 = 



f Q = 0 


( 53 ) 


We -will let the drag parameter CpA/m take on values of 1 and 10. The 
equations of motion are then solved for four cases: 


Case 1 = 0; = 0 

Case 2 \i / 0; = 0 

Case 3 i-L = 0;ft e /0 

Case 4 [i / 0; fte 4 0 


spherical gravitational field; 
nonrotating atmosphere 

nonspherical gravitational field; 
nonrotating atmosphere 

spherical gravitational field; 
rotating atmosphere 

nonspherical gravitational field; 
rotating atmosphere 


These cases permit study of the effects of oblateness and atmospheric 
rotation, at least for the above initial conditions. 


Accuracy of calculations.- The principal factors influencing the 
accuracy of the calculated results are the interval size and the number 
of significant figures carried in the calculations. Eight significant 
figures are carried throughout. The density was obtained by exponential 
interpolation in the ARDC density table and is not accurate to eight 
significant figures, although the density calculations repeat consistently 
to eight figures. With respect to interval size, the results of interval 
sizes in Ap of jt/8 and tc/ 32 are shown in the following list. This list 
gives certain quantities as calculated at cp = 13^ ^ or CpA/m = 1, case 1, 
initial altitude 120 miles : 



Acp = tc/ 8 radians 

Acp = rt/32 radians 

£, ft 2 /sec 

55. 084587x10 10 

55. 084366x10 10 

V 

0.21953070x10" 7 

0.23739280x10“ 7 

w 

- 0 . 12332460 X 10 " 4 

-0.12347577X10" 4 

t , sec 

34452 . 44o 

34452.uo 

r , ft 

21554337 

21554164 

r, ft/sec 

- 0.32462706 

-0.32502373 

A, radians 

4o. 840698 

4o. 840779 

E 

0 . 12702521 x 10" 4 

0.12717994x10" 4 

Vcp, ft/sec 

25556.150 

25556.253 

r, ft/sec 2 

-0.68629226xl0“ 6 

-0.71082937x10 " 6 

h, ft 

6 27909.00 

627736.OO 

P 

0.97086900 

0.97087600 
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The quantities are sufficiently close for most engineering purposes. 
Note that the altitude difference is 173 feet, which represents about a 
0 . 03 -percent error. Similar calculations were performed for case 2 with 
the same general results. On the basis of these results, it was decided 
to use tc/8 as the basic interval for the calculations. 

Periodic variations .- In figure 3( a ) are plotted the altitude and 
eccentricity variations for one period for CpA/m =1, cases 1 and 3* 

The loss in altitude with atmospheric rotation is slightly less than 
without rotation because the dynamic pressure and hence drag is decreased 
about 12 percent by rotation for eastward launchings. The path is a 
spiral for which the eccentricity, although small, has a strong periodic 
component. 


The principal effects of oblateness during a revolution are illus- 
trated by comparison of figures 3( a ) and 3( c )* The differences in the 
scales of the ordinates should be noted during the comparison. Thus the 
effect of atmospheric rotation is evident in figure 3( a ) hut not in 3(c). 
It is of interest that the altitude difference during one revolution is 
about 66,000 feet with oblateness, compared to about 800 feet without it. 
To explain this oblateness effect, we make use of the full gravitational 
field for an equatorial orbit as given by equation (10) 


G r 

m 


KM 



( 5*0 


The gravitational field is increased by 6p, or about 0.l6 percent. Since 
the satellite was initially in an equilibrium circular orbit with p = 0, 
suddenly "turning on" the oblateness at the initial altitude leaves the 
satellite with a velocity deficiency, just as if a tangential retrorocket 
had been fired. Since r Q = 0, the satellite location immediately turns 
into an apogee point; and, neglecting drag, the path becomes elliptical. 
From an energy consideration, neglecting drag, it is shown in appendix G 
that the altitude difference between perigee and apogee is 

r a " r p = 2r a e J € = ^ (55) 


We find that oblateness for equatorial orbits thus causes a difference 
in altitude between apogee and perigee of r a - rp = 66,700 feet. With- 
out oblateness but with drag figure 3( a ) shows an altitude loss for half 
a revolution of about k-00 feet. The difference in altitude between 
perigee and apogee is almost exclusively due to oblateness for the pres- 
ent value of CpA/m. In fact the calculated value from equation (55) of 
66,700 feet for oblateness alone is in very good accord with the value 
from the complete calculations of 67,081 feet for the combined effects of 
oblateness and drag. It is also clear that the approximate figure of 
67,000 feet is not sensitive to initial altitude since drag is not 
significant and the oblateness effect as calculated from equation (55) 
is insensitive to changes in altitude for near satellites 0 With regard 
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to the net loss of altitude per revolution, oh Lateness has the effect in 
the present case of doubling this loss. This specific result applies to 
the altitude range for which d(log p)/dh is the same as for the present 
example . 

For Cp A/m of 10, the altitude changes ani the eccentricity are about 
10 times those for CpA/m = 1, but there are ns qualitative differences. 

The effect of the rotation of the atmosphere has a small effect in all 
cases. 

Secular trends; limitations of calculation . - The manner in which the 
altitude and eccentricity vary from revolution to revolution, the secular 
trends, is illustrated in figure 4. For the range of figure 4(a) the 
decrease in altitude per revolution is a nearly linear function of the 
number of revolutions. However, for CpA/m of 10, figure 4(b), the rate 
of decrease of altitude becomes greater at lower altitudes since large 
density increases occur. While the altitude variations are as we would 
expect, the secular variations in eccentricity are of particular importance. 
For CpA/m of 1, figure 4(a) shows the eccentricity variation for 2 cycles, 
and then continues the envelope for 10 cycles, the limit of the calcula- 
tion. However, for CpA/m of 10, the eccentricity shows a divergence in 
the range of calculations, and the calculations in fact break down. 

To see how the calculations break down, let us examine the variation 
of eccentricity with cp as shown in figure 5* (For the equatorial orbit 
being considered cp and A differ only slightly.) Actually the oscilla- 
tory variation in E is no longer significant, as in figure 4(b), but E 
has what appears to be a nearly vertical tangent. The instantaneous 
ellipse is therefore undergoing very rapid changes in eccentricity. The 
interval size is too coarse to follow the rapidly changing curvature. By 
reducing the interval size, and increasing the number of figures carried 
in the calculation, the range of the calculatiDns can be increased. How- 
ever, this method is inherently unsuited to calculation of the terminal 
phase of the trajectory for the following reasDns: In the terminal phase 

of the trajectory the satellite descends verti 2 ally or nearly vertically 
along a radius vector. For such motion there is little or no change 
in cp, the independent variable. Therefore, a better-conditioned indepen- 
dent variable, such as time, should be used. 

Regression of line of nodes and movement pf line of apsides . - The 
solutions for the movement of the line of nodes and the motion of the 
line of apsides for the case of oblateness but no drag are known. The 
variation in these quantities per revolution cm be obtained by integrat- 
ing the equations for d0/dcp and db/dcp from p = 0 to 2jt on the basis 
that departures from the basic ellipse are small. The following results 
are taken from reference 13. 


1 ?_h 
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(56) 


where V c is the circular orbital speed at radius r for no oblateness 
(eq. (53))- It is to be noted that the first equation, strictly speaking, 
applies to the secular change for one nodal period , 4 that is, for the time 
to go from one ascending node to the next. It happens that drag forces, 
being in the plane of the orbit, have only a very small influence on the 
regression of the line of nodes, as shown by the following list. 



0 , radians 

2* 

°D a , 0 

= 1: case 2 

m * 

CpA 

- 10; case 2 

m y 

Equation (56) 

1 

-0.009689 

-0.009695 

-0.009696 

2 

-.01938 

-.01941 

-.01939 

3 

-.02907 

-.02914 

-.02909 

4 

-.03876 

-.03892 

-.03878 


The movement of the line of apsides (the major axis of the instantane- 
ous ellipse) is vitally influenced by the drag* It is first desirable to 
note what happens in the absence of drag for an equatorial orbit. The line 
of nodes, which can be visualized only for orbital planes away from the 
equatorial plane, moves backward against the motion of the satellite at 
0.555° P er revolution. The corresponding rate for the moon is about 
1.5° per revolution. The line of apsides moves forward in the direction 
of the satellite at twice this rate. As a result, the line of apsides 
moves around the equator with respect to an inertial framework at precisely 
the same rate that the line of nodes moves backward. 

With the introduction of drag the line of apsides tends to move 
around the orbit with variable lag, at the average speed of the satellite. 
To show clearly the motion of the line of apsides consider the angle 5 
given by equation (48) 

& = tan" 1 £ (57) 

The angle 5 gives the position of perigee 5 measured from the line of 
nodes. A plot of w versus v shown in figure 6(a) illustrates how & 
starts out at ^ for cp = 0 for the case shown and increases steadily 
thereafter. If 5 leads or lags cp by a constant amount at all times, 

4 The nodal period is sometimes referred to as the synodic period. 

5 The position of perigee is sometimes referred to as the minor apsis. 
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then a curve of (5 -cp) versus cp would, he a straight line. The parameter 
(6 -cp)/2 it is shown versus cp/2jt in figure 6(h) for the present case. 

The variable lead of 5 over cp approaches ji/ 2 after several revolu- 
tions. The principal conclusion to he drawn from these results is that 
the very small motions of the line of apsides due to ohlateness alone are 
completely masked when drag is included because the line of apsides moves 
with average speed of the satellite. The concept of the line of apsides 
becomes useful now as a measure of lead or lag of the angular position of 
the major axis of the instantaneous ellipse from the angular position of 
the satellite. 


Second Illustrative Example 


As a second illustrative example, consider a satellite launched from 
the equator in an orbital plane inclined at 65° to the equatorial plane. 
The initial height above the equator is taken as 80 miles, and the satel- 
lite is launched at circular orbital speed. Tie initial altitude was 
chosen low enough to obtain impact in about a revolution or less. The 
launching velocity was not always horizontal as in the previous example, 
but radial velocities were introduced as follows: r 0 /(Vq>) 0 = 0, -0.01, 

- 0 . 05 , - 0 . 10 . Values of CpA/m of 1 and 10 are included. The equations 
of motion were integrated numerically for the four cases considered in 
the first illustrative example. An attempt was made to carry all cases 
to impact. Impact altitude is taken to correspond to 1000 feet altitude. 
The calculations have been made using geographic coordinates, but were 
checked against calculations in one case using the elliptic orbital 
elements method. 

Accuracy of calculations .- The accuracy oi the calculated results in 
this example can be assessed by several means. First, we can compare 
calculations on the basis of geographic coordirates with those on the 
basis of elliptic orbital elements 0 We can also vary the time interval 
used in the calculations based on geographic coordinates. Pursuing the 
first comparison, consider a satellite with CpA/m = 1, r 0 /(^cp) 0 = 
case 1, for cp = 6.037 radians. The orbital quantities for these condi- 
tions are compared in the following list. 



31 


Orbital elements, Geographic coordinates, 
Acp = tt /64 radians At = 30 seconds 


cp, radians 

6.0377500 

6.0377582 

V 

- 0.003810 776^ 

-0.0038172179 

w 

0.0030982814 

0.0031003921 

t , sec 

5010.6056 

5010.6056 

r, ft 

2. 1258168x10 7 

2.1258139XIO 7 

r, ft/sec 

-54.238659 

-54.265940 

A, radians 

6.1777185 

6.1777180 

E 

0 . 0049667594 

0.0049733320 

Vq>, ft/sec 

2.5675599X10 4 

2.5675529XIO 4 

r, ft/sec 2 

-0.13947319 


h, ft 

0.335174x10 s 

0.335146x10 s 

P 


0.98883584 

^ , radians 

-0.22203323 

-0.22202809 

a, radians 

1.1344640 

1.1344641 


It is noted that the two sets of calculations are in good agreement, 
indicating that no gross mistakes have occurred either in analysis or 
calculation. The interval size for the orbital elements is the 128 th 
part of a revolution, while that for the geographic coordinates is 
essentially the l80th part of a revolution* 

The foregoing comparison is made for a satellite descending from 
80 miles altitude to about 63 miles altitude in one revolution. When 
the calculations were performed with the 30-second interval, it became 
clear on numerical grounds after some point in time that the interval 
size was too large. All ealculative cases using a 30-second interval 
exhibited this behavior before impact. Thus while an interval size of 
30 seconds is satisfactory for the initial phase of the trajectory, it 
is not adequate for the terminal phase. Some time before the calcula- 
tions exhibited inaccuracies with the 30-second interval, the interval was 
switched to 3 seconds. The calculations were then continued to impact. 

Let us now compare the trajectories calculated using the foregoing 
method based on two interval sizes with the trajectory calculated using 
a 3 “Second interval size all the way. Since we are principally interested 
in the point of impact, the following comparison yields a good idea of 
how accurate the point of impact is. The example considered for this 
comparison corresponds to CpA/m = 1 , case 1 , r 0 /( v cp ) 0 = the left- 

hand column a 3~ sec ond interval was used to compute the trajectory all 
the way to the point of impact. In the right-hand column, the 30 -second 
interval was used to t = 5760 seconds and the 3-second interval from 
then until impact at 5991 seconds. 
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Interval size, sec 

3 only 

30 and 3 

t , sec 

5991 

5991 

cos cp 

0.86217703 

0.86456735 

r, ft/sec 

166.43243 

167.09205 

A, radians 

6.5265953 

6.5241574 

v cp> ft/sec 

0 

0 

h, ft 

327 

592 

i|r, radians 

0.47702802 

0.47302309 


In this case the satellite travels slightly more than a complete revolu- 
tion around the world. The difference in latitude at the end of 5991 sec- 
onds is 15.87 miles and the difference in longitude is 8.58 miles. Since 
the range is over 25,000 miles , these accuracies for the impact point 
were considered satisfactory. Hence most of the calculations were made 
on the basis of a 30 -second interval followed by a 3 -second interval in 
the t e rmi nal pha s e • 

General features of the trajectories .- To show the general features 
of the trajectories , consider figures "7~(a) , 7 (t), and. 7 ( c ) constructed 
for Cjyk/m = 1, r 0 /(Vcp ) 0 = 0. The altitude variations with time shown 
in figure 7 (a) exhibit several important effects. The general waviness 

is due to the nonspherical figure of the earth. At t = 0 the satellite 

starts at the equator, and at about 1200 seconds the satellite reaches 

its maximum north latitude for which the earth radius is least for the 

orbit and the altitude is correspondingly greater. Subsequent passes 
over the equator through maximum and minim-urn latitude cause further 
bumps. In the cases including the gravitational effects of oblateness, 
the time of flight is significantly shorter than without oblateness 
effects. The basic reason for this behavior has already been discussed 
in connection with the equatorial orbits. Oblateness causes the satel- 
lite to descend lower into the atmosphere. The resulting higher drag 
thus reduces the flight time as shown. 

The latitude -longitude variations of the satellite for the four cases 
are given in figure 7(t). The paths for the four cases are essentially 
the same until the time atmospheric drag Initiates entry. At this time 
the latitude and longitude are essentially frozen. What this means . 
generally is that at a given time the satellite for the four different 
cases has nearly the latitude and longitude given by the Kepler ian solu- 
tion, but the altitudes differ signif icantly . As a result, for the four 
different cases the satellite enters the final constant longitude -latitude 
phase at different times. 

During the first part of the trajectory the satellite not only 
follows the same latitude -longitude path for all four cases, but it also 
appears at a given longitude and latitude at the same time. However, 
near the very end of the trajectory the satellite decelerates rapidly 
just prior to turning down into the atmosphere. During this phase the 
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satellite appears at a given latitude and longitude somewhat later than 
.for those cases for which the satellite has a greater altitude. Included 
on the curves of figure 7(a) are ticks which indicate where the eccen- 
tricity nearly reaches unity. These ticks correlate closely with the 
position where the latitude and longitude become constant. 

The eccentricity variations of figure 7(c) shed some light on the 
motion. The eccentricity started at zero because of the particular 
initial conditions taken in the present calculation; namely, those cor- 
responding to a circular satellite orbit (see eq. (Bll)). The eccen- 
tricity remains small for the first part of the trajectory and shows 
characteristic waviness because of the nonspherical figure of the earth. 
During the final half of the trajectory the eccentricity rapidly rises 
nearly to unity and remains there during the terminal phase. The value 
of unity is associated with a vertical path as shown by equation (B7)* 

For cases 1 and 2 the eccentricity remains essentially unity until impact, 
in most instances. In other instances, however, the machine calculation 
of eccentricity becomes erratic. Actually , for a vertical path the 
eccentricity is a derived quantity which has no particular significance, 
so that the erratic behavior of E does not reflect on the accuracy of 
the trajectory calculations. 

Impact point . - The calculations of this illustrative example also 
shed some light on how the oblateness component of gravitational field 
influences the impact point, as well as how atmospheric rotation influences 
it. Let A x ; , A 2 ; \jr 3 , A 3 ; A 4 be the impact coordinates for 

cases 1, 2, and. > respectively. Then the effects of oblateness on 
the impact point with no atmospheric rotation are >|r 2 “ A 2 - A x and 
with atmospheric rotation are i|r 4 - 4r 3 , A 4 - A 3 . These quantities are 
tabulated in table I as a decimal part of the total range from the 
assumed initial point to impact. The coordinates A 4 , i|r 4 of the impact 
point are also tabulated together with the great circle distance, s, 
between 0, 0 and A 4 , ty 4 . These latter quantities as calculated in 
radians were multiplied by a, the equatorial radius, to convert to 
miles. The ranges are as measured in the inertial system XYZ and are 
not those for an observer on the earth. The earth ranges can be deter- 
mined with the help of the tabulated flight times. The first important 
point is that oblateness causes large percentage errors in range for 
re-entry at zero or very small angles. These errors are simply the 
calculated differences in impact point between those for p equal to 
zero and for p not equal to zero. As previously mentioned, the initial 
velocity is too low to launch the satellite into a circular orbit when p 
is not zero. 

In the foregoing discussion the effects of oblateness have been 
taken as differences in the calculated results due to neglecting the 
oblateness gravitational terms in the equations of motion for the same 
initial conditions. It appears desirable to take some account of oblate- 
ness in the initial conditions. For equatorial orbits this is easily 
accomplished by making V<p equal to that for a circular orbit. Without 
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oblateness the value of V™ to achieve a circular orbit is V c , but 
with oblateness the initial Vq> should be 


Vcp = V co = 



Let us now differentiate the following three examples: 
Example Initial V f 


*cp 


p m 

equations of motion 


A 

B 

C 


Vcp = V, 
Vcp = V 


Vcp = V co 


Not included 
Included 
Included 


For these three examples trajectories were computed to impact for the 
following initial conditions: 


a r 


= 0 " 


j h-o 


= 80 miles; C]>A./m = 1 ; 

o 


r o ” 0 j 9 q ~ ^o - to ~ ^ 
The following results were obtained: 

Example A Example B 


Example C 


h, ft 

151 

270 

193 

t , sec 

4,058 

2,742 

4,017 

Vcp, ft/sec 

lo " 11 

10 _1G 

10 " 10 

r, ft/sec 

165.883 

166.310 

166.121 

A, radians 

4.4820 

2.9011 

4.4338 

A or s, miles 

17,750 

CO 

o^ 

17,570 

The foregoing results 

are as one 

might anticipate 

. For ’ 


^cp - V c 

the range is reduced 6,252 miles by including oblateness in the equations 


of motion, but with V, 


<P 


V, 


co 


the range is reduced only l 8 o miles. 


If the inclination of the orbital plane to the equator is changed 
from zero, a purely circular orbit with oblateness is not possible. As 
a matter of curiosity the value of Vcp for a nonequatorial orbit was 
changed from V c to V co to see how the range was affected. The case 
investigated corresponds to Oq = 65 °, Cjyk/m = 1 , r 0 = 0 , and h Q = 80 miles. 
The ranges are tabulated for the same conditions as examples A, B, and C. 
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Example A; s = 26,987 miles 
Example B; s = 18,73 8 miles 
Example C; s - 38,126 miles 

The ranges are great circle ranges. It is noted that including p 
in the equations of motion but not adjusting the initial value of V<p 
reduces the range by about 30 percent. If the initial value of is 

boosted to V co , the range is greatly increased for the present orbit 
even though p is included in the equations of motion. 

It is clear that for the present initial conditions increasing 
CpA/m from 1 to 10 diminishes the effect of oblateness on impact point. 
This result is in accordance with the fact that increasing CpA/m 
decreases the importance of gravitational forces compared to drag forces. 
Increasing the initial downward velocity also decreases the proportionate 
error in impact point. This effect is similar to decrease in miss dis- 
tance due to errors in launching speed as the downward launch angle is 
increased. 

For manned re-entry for which CpA/m is of the order unity and the 
entry angle is of the order of a degree or less (r/Vq) < - 0 . 02 ) to limit 
normal accleration, oblateness has an influence of several hundred miles 
in latitude and longitude on the point of impact. Including or neglect- 
ing atmospheric rotation does not significantly influence this result. 

The influence of atmospheric rotation on the impact point is gener- 
ally not so large as that of oblateness. The influence of atmospheric 
rotation is represented by ^ 3 -^, A 3 - A x , or \[r 4 -^ 2 , A 4 - A 2 . These 
quantities are tabulated in table II in the same manner as table I. 

As the downward launch angle is increased, the percentage effect of 
atmospheric rotation on range generally increases. Such an effect is 
due to the fact that the satellite spends more of its time in the lower 
atmosphere where the air density is higher. Also, as CpA/m is increased, 
the drag due to atmospheric rotation assumes more importance and increases 
the effect of atmospheric rotation on impact point. As the satellite 
approaches its impact point, its vertical velocity usually is small 
compared to the rotational speed of the atmosphere. The atmosphere thus 
drags the satellite around with it at constant ^ and increasing A. As 
seen by the observer on the earth, the satellite would descend vertically 
except for a small slippage between the satellite and the atmosphere. 

For a value of CpA/m of 1 and for £ 0 /(Vq )) 0 of -0.01, such as might 
be used for manned entry, a 150 -mile change in impact point is due to 
atmospheric rotation for the present initial conditions. 
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RELATIVE IMPORTANCE OF VARIOUS TERMS 
IN EQUATIONS OF MOTION 


From the complete trajectories presented Ln the previous calculated 
examples, it is possible to examine the order of magnitude of the various 
terms occurring in the equations of radial and tangential motion. The 
resulting information will show what physical terms are important for 
various phases of the trajectory and, hence, will suggest simplifications 
permissible for analytical work in particular. Also, the information will 
show where in the trajectory the assumption stated in equation ( 3 ) becomes 
valid; that is, where the entry or terminal phase begins and hence where 
the solution of Chapman (ref. k) is valid. 

The parameter which has particular significance for dividing the 
trajectory into various phases is the ratio of tangential velocity Vq> 
to circular satellite velocity V c . To obtain a number whose logarithm 
is not -00 when Vq> = 0, let us use 1 - Vq)/V c as the parameter. In fig- 
ure 8(a) the variations with 1 - Vcp/Vc are shown of the terms in the 
tangential equation of motion 


• i’Vq) % Gcp 

V -f - — + — 

t r mm 


(58) 


The term rVcp/r is less than 10 percent of Vrp when Vcp is less than 
about 99’percent V c . This result gives a quantitative measure of when 
the assumption of equation ( 3 ) is met. Figure 8(b) was constructed to 
show that the figure of 99 percent does not change when the drag parameter 
CpA/m is increased from unity to 10. 

Figure 8(b) exhibits phenomena not manifest in the range of fig- 
ure 8(a). First it is seen that Keplerian motion characterised by 
Vcp = -rVcp/r is never realized for this example. The sequence of events 
is interesting to examine. At time zero the drag is in equilibrium with 
the tangential acceleration force and the satellite slows down. However, 
although drag initially makes Vcp negative, tie satellite tends to speed 
up as it drops in altitude until Vcp is zero. At this point the drag is 
in equilibrium with the acceleration force due to rVcp/r. The satellite 
speeds up as it drops in altitude until a maxi mum value of Vcp is 
reached, when Vcp is again zero. Thereafter the satellite decelerates 
tangentially at an increasing rate, and the motion is in accord with the 
solution of Chapman. For CpA/m = 1 figure 8(:) shows that setting 
r D = -0.l(Vcp) o does not alter the range of applicability of Chapman* s 
assumption. 


With regard to the equation of radial motion 



m 


(59) 
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the magnitudes 
unity and r Q 


of the four terms are shown in figure 9 ( a ) for CjyA/m of 
of zero. It is interesting to note that the gravitational 
and centrifugal forces are in balance until Vcp is about 0.9 V c . There 
is then an interval where four terms are important. For V<p near zero, 
the drag and gravitational terms are nearly in balance. However, the 
radial acceleration term r cannot be ignored If accurate terminal 
velocities are desired, since it is about 30 percent of the other two 
terms. If CpA/m is increased to 10 from unity (fig. 9(^0 )> or if the 
initial launch angle is 0.1 radian downward (fig. 9(c)), the radial 
acceleration still cannot be neglected except for rough calculations. 

A patching technique is useful to establish a rough complete trajec- 
tory. For Vcp above about 0.99 V c the trajectory can be approximated 
by the method of appendix A. For Yep equal to about 0.99 Vc the solu- 
tion can be joined to that of Chapman, as described in reference 4, and 
continued down to the point where the flight path is nearly vertical. 
Although the equation of Chapman is still valid for vertical flight, the 
numerical solution of his equation loses accuracy because of a singularity. 
A solution for the vertical part of the trajectory is given in 
reference l4. 

Let us examine the contribution of atmospheric rotation to the drag 
term since this has important implications concerning the adequacy of 
approximate two-dimensional theories which usually neglect atmospheric 
rotation. This question assumes significance in the terminal phase of 
the trajectory. Let us examine the flight path angles as seen by an 
observer in an inertial framework and by an observer on the earth. The 
flight path angle in the inertial framework is simply given by 


tan n = v" 

v cp 

For an earth observer, r is unchanged but Vcp is decreased for eastward 
motion by the component of the earth r s rotational speed in the orbital 
plane. This component of speed is r cos ^O e sin p or rft e cos a. The 
flight path angle, / e , as seen by an earth observer is simply given by 


tan 7 e = 


Vcp - rftecos a 


To illustrate the influence of earth rotation and observational position 
on flight path angle, figure 10 has been prepared. The initial conditions 
are taken to be 

a 0 = 65°; (Vcp) 0 = (V C ) Q ; r Q = 0; h Q = 80 miles; 

CjyA/m =1; d 0 = A 0 = \|/ Q = 0° 


The first curve to which attention is called is the plot of 7 ^ for 
case 1, \i = 0, and Q, e = 0. As the satellite approaches impact for this 
case, the flight path approaches a nearly vertical condition as seen by 
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an observer in the inertial framework. If we now include atmospheric 
rotation, (and oblateness), the variation of for case 4 is obtained 

as shown. The behavior of 71 in this case is distinctly different since 
it now reaches a maximum near 30° and rapidly decreases thereafter. What 
has occurred in the terminal phase is that the radial speed has been 
reduced to several hundred miles per hour by air drag, but the horizontal 
speed Vcp has been increased to a much greater value by atmospheric rota- 
tion, The atmospheric rotation causes the trajectory to curve almost 
directly eastward at constant latitude \|r with a horizontal speed equal 
very nearly to rQ e cos \|r. At the equator this speed is about 6 percent 
circular orbital speed for near earth satellites. Thus, approaching 
impact, the satellite has a nearly constant value of V^/Vq less than O.06. 

The distinctly different characteristics cf the trajectory in the 
terminal phase with and without atmospheric rotation raise the important 
point whether two-dimensional theories neglecting atmospheric rotation 
are really applicable to the terminal phase. The question can be 
answered in the affirmative, provided the correct interpretation is given 
to the two-dimensional theories. In this connection the flight path 
angle as seen by an earth observer has been plotted in figure 10 for 
case 4 which includes oblateness and rotation. Not unexpectedly, it 
turns out to be in close accord with for case 1, no oblateness or 

rotation. The following interpretation is given to this result: Two- 

dimensional theories neglecting atmospheric rotation (but including drag) 
yield trajectories which tend to be nearly vertical in the terminal phase 
provided CpA/m is not too small compared to unity. These theories can 
be applied with good accuracy to three-dimensional trajectories including 
atmospheric rotation if the results are interpreted to apply to the 
motion as seen by an earth observer. In fact, vith this interpretation 
it would be a mistake to include atmospheric rotation in the two- 
dimensional theory. 


COMPARISON OF EQUATORIAL TRAJECTORY WITH TRAJECTORIES 
OF APPROXIMATE TWO-DIMENSIONAL THEORIES 


It is of interest to compare approximate t vo-dimensional analytical 
results with an equatorial trajectory as calculated numerically by the 
present method. For this reason a special equatorial trajectory was 
computed neglecting atmospheric rotation and oblateness effects, factors 
not usually considered in two-dimensional theories. The initial condi- 
tions for the trajectory are 

h Q = 80 miles; (V (p ) Q = (V c ) Q J r G = 0; 

= ~ = G-O = 0 

The trajectory was calculated for a 1-second interval using geographic 
coordinates. The values of various quantities luring re-entry are 
listed in table III. 
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The approximate theories compared with the present method include 
that of appendix A, the method of Chapman (ref. 4), and an extension of 
the results of Linnell (ref. l4). These theories are applied, respec- 
tively, to the initial part of the flight, the middle part of the flight, 
and the terminal phase of the flight. 

Let us first consider that part of the flight from Vqj/Vc from 
about 1.000 to about 0.995^ during which most of the range occurs. The 
small variations in Vm and r permit the simple approximate solution of 
appendix A. The variations of h and r with t obtained by this method 
are compared in figure 11 with those given in the foregoing table. Solu- 
tions are given for two time intervals in the calculation, a short interval 
and a long interval as described in appendix A. For the short time 
interval the values of r and h are in fair accord with the tabulated 
values out to larger values of t than for the longer time interval. 

Consider now the phase of the trajectory with Vqj less than 
99-5- percent circular orbital speed. This is the region where the 
rVq)/r force is negligible and the solution of Chapman applies. This 
solution is expressed in the form of a parameter Z tabulated as a 
function of u 



The Chapman solution considered here is that for u = 0.995 and 7 = -0.5°, 
and is the available one most nearly approximating the present calcula- 
tions. For the initial values of u and 7, the solution of Chapman gives 
a value for Z. From this initial value of Z and the value of CpA/m 
of unity, the initial value of the density can be calculated from equa- 
tion ( 60 ). The values of h, r, Vm, and A are compared in figure 12(a) 
and the value of 7 is compared in figure 12(b). The initial altitude 
obtained from the Chapman solution is slightly less than that obtained 
from the present solution. Some of this difference is due to the slight 
differences in 7 exhibited in figure 12(b). Generally speaking, the 
solutions are in good accord. The tendency of the solutions for 7, Vm, 
and h to interlace is probably due to slight differences in the atmos- 
pheric altitude-density relationships assumed in the two methods. 

Chapman has used an exponential atmosphere in his work, while the present 
work is based on the ARDC atmosphere. 

In the terminal phase, the present tabulated solutions of Chapman do 
not continue entirely to impact but stop at a value of u = 0.025. For 
u = 0, vertical flight, the solution of Linnell, reference l4, is 





available. In an unpublished investigation Elliott D. Katzen of Ames 
Laboratory has adapted the solution of Linnell to slight departures from 
vertical flight. His calculated values for u less than 0.025 agree 
well with the present values. 


CONCLUDING REMARKS 


The principal purpose of this paper has been accomplished; namely, 
to present equations for calculating complete trajectories of earth satel- 
lites from outer space to the ground under the influence of air drag and 
gravity, including oblateness 6 and to apply these to several cases of 
entry from circular orbits. 

Equations of motion based on an "instantaneous ellipse" technique, 
with polar angle as the independent variable, were found suitable for 
automatic computation of orbits in which the trajectory consists of a 
number of revolutions. This method is suitable as long as the satellite 
does not enter the terminal phase. In the terminal phase of the trajec- 
tory, equations of motion in spherical polar coordinates with time as the 
independent variable were found to be suitable. 

In the first illustrative example, the effects of the oblateness 
component of the earth's gravitational field ard of atmospheric rotation 
were studied for equatorial orbits. The satellites were launched into 
circular orbits at a height of 120 miles, an altitude sufficiently high 
that a number of revolutions could be studied. The importance of the 
oblateness component of the earth's gravitational field is shown by the 
fact that a satellite launched at circular orbital speed, neglecting 
oblateness, has a perigee some 67,000 feet lower when oblateness forces 
are included in the equations of motion than when they are not included. 
Also the loss in altitude per revolution is double that of a satellite 
following an orbit not subject to oblateness. The effect of atmospheric 
rotation on the loss of altitude per revolutior was small. As might be 
surmised, the regression of the line of nodes as predicted by celestial 
mechanics, equation (56), is unchanged when drag is included. It is 
clear that the inclination of the orbital plane to the equator will be 
relatively unaffected by drag for no atmospheric rotation since the drag 
lies in the orbital plane in this case. With "the inclusion of atmospheric 
rotation it was found that the inclination of the plane changed about 
10 -6 radians per revolution. Thus the prediction of the position of the 
orbital plane of an earth satellite is not complicated by the introduction 
of drag. The line of apsides, which without drag but with oblateness 
moves slowly in space, tends to move with the i atellite when drag is 
included in the calculations. As a result the usual linearized solutions 
based on oblateness alone must be basically altered when drag is included 
to take into account the rapid movement of the line of apsides. 

6 No attempt was made herein to take into account gravitational 
anomalies or surface cross winds . 
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In the second illustrative example, the final revolution was 
calculated to impact for a number of trajectories in an orbital plane 
inclined at 65 ° to the equator. Of particular interest are the large 
effects that the oblateness gravitational field and atmospheric rotation 
have on the impact point. For a value of Cj)A/m of unity, and for an 
initial downward angle at 80 -miles altitude of 0.01 radian, such as might 
be utilized for manned re-entry, oblateness had an influence of about 
300 miles in the impact point, and atmospheric rotation had about a 
150 -mile influence . 

From the complete trajectories calculated automatically, it was 
possible to examine the relative importance of the various terms in the 
equations of motion. For the equation of tangential motion, the term 
proportional to rVq)/r is negligible, as long as the local value of Jcp 
is less than 99 percent of V c . This result indicates that for less 

than 0.99 V c the equation of Chapman (ref. 4) can be used. For the equa- 
tion of radial motion, the radial component of gravity and the centrifugal 
force dominate the motion for values of Vq> near V c . However, for 
Vcp « V c > where the trajectory is nearly vertical, the radial component 
of gravity and the drag dominate the radial motion, but the radial 
acceleration is not generally negligible. 

It was found that two -dimensional solutions neglecting atmospheric 
rotation can be used to approximate three-dimensional solutions with 
atmospheric rotation. In this connection, two-dimensional theories must 
be interpreted as being viewed by an observer on a rotating earth. 

Several gaps exist in the solutions available for studying the 
dynamics of earth satellites. First, to the authors 1 knowledge, no 
linearized theory exists for predicting the periodic variation of the 
elliptic elements during circularization of the ellipse or during spiral 
decay, taking into account drag. This linearized solution would have 
to take into account the fact that the line of apsides tends to move with 
the satellite. Second, a missing ingredient in the accurate automatic 
computation of impact points is precise knowledge of the variation in 
atmospheric density with latitude, season, or time of day. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., July 15, 1958 
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APPENDIX A 

APPROXIMATE SOLUTION FOR TWO-DIMENSIONAL TRAJECTORIES NEAR 
CIRCULAR ORBITAL SPEED TOGETHER WITH NUMERICAL EXAMPLE 


Since a considerable part of the range of a satellite re-entering 
the atmosphere is traversed at a tangential speed V<p nearly equal to 
circular orbital speed, it would be helpful to have an approximate solu- 
tion for this case. In what follows we assume that the percentage 
changes in Vq> and r are small, but the perce itage change in air 
density p may be large. We will find an approximate solution r = r(t) 
subject to the initial conditions r = r Q , r = r 0 at t = 0. (The initial 
value of *r 0 can be .calculated from ro and r c*) 


Consider the equations of motion in the following form: 


_ V = 

r m m 


• . Vm Dm Gm 

v c? + r — = — + — 


(Al) 


If we neglect oblateness, then Gy Is zero. Also, the flight path angle 
is nearly zero so that Dr can be neglected as shown in figure 9- The 
equations of motion with no atmospheric rotation are 


r = 


*cp 


KM 


Vcpr + rV 9 = - 2 P V 



> 


(A2) 


Since the only variable in the tangential equation of motion with any 
appreciable percentage change is p, we can integrate the equation as 
follows 


( r V = < r V 0 • I 


(A3) 


The quantity p is the mean density during the time interval between 
t = 0 and t = t, and (rVcp) 0 is evaluated at t = 0. Equation (A3) gives 
the time range of change of area sweeping. Since Vy and r are slowly 
varying functions of time, the product rVy aLso is. 
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To obtain the altitude as a function of time necessitates an 
integration of the radial equation of motion 

•• = _ KM (rVcp) 2 
r 2 + r 3 

Equation (A3) together with the relationship 

r = r D + r Q t + 0(t 2 ) 
put equation (A4) into the form 

r = r Q - kt + 0(t 2 ) 


(A4) 


(A5) 
(A 6) 


wherein 



Integration of equation (A6) yields 



3r 0 + 



r = r Q + r 0 t - ^|- + 0(t 3 ) 
r = r o +i’o t +^o \ + °(t 4 ) 


(A7) 


(A8) 


To obtain the altitude time curve for a particular case, we first 
choose a value of p* for the first step and evaluate k, The value of 
r versus t can then be established from equation (A8). The curve of 
r versus t also establishes a curve of p versus t for a given atmos- 
phere, The time t x for which p - satisfies the following relationship 
can easily be found 

P = f [ X p dt 

z xJ 0 

Since p is the average value of p between t = 0 and t = t x , the 
values of T ly r x , and r x calculated from equations (A6) and (A7) for 
t 1 should be accurate. The values of r x , r-L, and r x are now used as 
initial values in the next step of the calculation. The next step is 
started by moving the zero of the time scale to t x , choosing a new value 
of p3 and calculating a new value of k. Thus the process is continued 
step by step. 


(A9) 
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As a numerical example of this calculation consider a satellite 
launched in an equatorial orbit in an easternly direction with an initial 
altitude of 80 miles. The initial conditions necessary to proceed with 
the calculations are 


t G = 0 sec 

r Q = a + h Q = 21,348,828 ft 
fo = 0 ft/sec 
(Vcp) 0 = 25,678.862 ft/sec 


CpA 

m 


1.0 


J 


(A10) 


At an altitude of 80 miles the density from reference 7 is p = 2.96x10 11 
slugs per cubic foot. For the first step in the calculation of the 
trajectory let us choose a value of p slightly greater than the initial 
density, as follows, 


p = 3*019x10 11 slugs/ft 3 
and evaluate equation (A7) for k. 


(All) 


k = (3. 019x10" 11 ) 


(25,678.862; 3 
(21,348,828) 


( 1 . 0 ) +0 


= 2.3945X10" 5 (A12) 

The considerations in selecting p* are discussed subsequently. 

With the value of k determined, equaticn (A8) is now used to 
establish a curve of r versus t. For t eqval to 0, 400, 800, and 
1200 seconds the respective values of r are found to be 


t = 0 sec r = 21,348,828 - 10 — (0]_ 

6 

= 21,348,828 ft 
t = 400 sec r = 21,348,573 ft 

t = 800 sec r = 21,346,785 ft 

= 21,341,932 ft 


t = 1200 sec 


r 
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These points also establish a curve of p versus t by the use of refer- 
ence 9 to obtain p for a given altitude above the earth f s surface. The 
corresponding values are 


t = 0 sec 

P = 2.96x10 11 slugs/ft 3 

t = 400 sec 

p = 2. 98xlO~ 11 slugs/ft 3 

t = 800 sec 

p = 3. 18x10' 11 slugs/ft 3 

t = 1200 sec 

p = 3. 78x10" 11 slugs/ft 3 


From this curve of p versus t the 
value of t x which satisfies equa- 
tion (A 9 ) can be easily found, as 
is shown in the accompanying sketch, 
by making the two shaded areas 
under the curve equal. Performing 
this integration for this numerical p 
example one finds that 

t 1 = 821 sec (A 13 ) 

Now, using the values from 
equations (A 10 ), (A12), and (A 13 ) 
and inserting them in equa- 
tions (A6) and (A8) yields the 
quantities 



♦ 

Sketch (c) 


r i 


21,348,828 


( 2 . 394 ^xl 0 ~ 5 )( 82 i) 3 

6.0 


= 21,346,619 ft 

(2.3949x10~ 5 ) (821) 2 


= -8.06996 ft/ sec 
r 1 = -( 2 . 3945 xlo" 5 )( 82 l) 
= -O.OI966 ft/sec 2 


Therefore, at t x = 821 seconds the quantities specifying the satellite’s 
position and motion are 
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n = 21,346,619 ft 
i*i = -8.06996 ft/sec 
= -O.OI966 ft/sec 2 
(Vcp)i = 25,678.862 ft/sec 


(A14) 


These values are now used as initial values in the next step of the 
calculation. The zero of the time scale is shifted to correspond to t x , 
and a new value of p must be chosen. The value of r x corresponds to 
an altitude above the earth's surface of 420,191 feet for which the 
density (ref. 7 ) is 3-205x10 11 slugs per cubic foot. As the average 
density for the next step, let us, therefore, choose a new value of 


p = 4. 30x10" 11 slugs/ft 3 (A15) 

With this value and equations (Al4) the calculation of k (eq. (A7)) is 
repeated and is found to be 


k = 2. 2452x10 " 5 


Now the curve of r versus t is determined from equation (A8) and the 
corresponding densities are determined from reference 9- These are 


), sec 

r, ft 

p, slugs/ft 3 

0 

21 , 346,619 

3 . 205 x 10" X1 

400 

21 , 341,579 

3 . 850 x 10" 11 

800 

21 , 331,956 

5. 390x10" 11 

1200 

21 , 316,314 

9 . 650 x 10 " 11 


Integrating graphically as before 


t 2 - t x = 959 sec 


and equations 


(A 6) and (A8) thus give 


r 2 = 21,326,540 ft 
r 2 = “37-24824 ft/sec 
r 2 = -0 ,04119 ft/sec 2 




(A16) 


(A17) 


Therefore, at t 2 = 1780 seconds, these quantities plus Vcp, obtained 
from equation (A3) specify the satellite's position and motion. 
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These quantities thus become the initial conditions for the next 
step and the calculation is repeated. 

Let us now consider the factors entering the choice of p. For 
entry from circular orbits , as in the present example , the value of p> 
should be 1 or 2 percent greater than the initial density, p Q . In the 
calculation of figure 11 for the long time interval, it was decided 
a priori to go to 3000 seconds in about four steps. The values of 'p/p 
for each step are shown in the following table: 


Time interval 

P/P 0 

Pi/Po 

0 - 821 sec 

1.020 

I.O83 

821 - 1780 sec 

1-342 

2.031 

1780 - 2425 sec 

I.966 

3.932 

2425 " 2875 sec 

3.125 

9.415 


Because the density changes slowly with altitude at first, a large initial 
value of p’/pQ would give a very large time interval. This is to be 
avoided since the present method is based on power series in time. Once 
the curve of p versus t is established in the first interval, the values 

of p/po for subsequent intervals to obtain given intervals in time can 

be estimated by extrapolating the curve. 

To study the effects of time interval, the calculation was also made 
in about eight steps instead of four. As expected, the calculation with 
more steps remains closer to the machine solution at large values of the 
time. In any particular case, it is best to do the calculation with two 

different time intervals to be sure of the range of accuracy of the 

calculations . 

It might be mentioned, in conclusion, that the present method can 
be applied to an atmosphere of arbitrary density-altitude relationship. 
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APPENDIX B 

APPROXIMATE EXPRESSION FOR ECCENTRICITY FOR 
VALUE OF Vqj SMALL COMPARED TO V c 


It is possible on the basis of energy considerations to obtain a 
simple expression for the eccentricity for values of Vcp/V c small 
compared to unity. The starting point is equation (4l) which relates 
the eccentricity to certain quantities as follows 


E 2 = 

pl 


From the definitions of L 0 and p we find 


i£ = lof ( ^ A\ ; 

pi KM \£ 0 2 / \IS KM 


(Bl) 


(B2) 


The length of the semimajor axis Z is related to the total energy for 
a circular force field. The kinetic energy for unit mass is V 2 /2 and 
the potential energy per unit mass is -KM/r. Thus the total energy E 
is 



(B3) 


Since the total energy depends only on the length of the major axis 
independent of the eccentricity, we have (ref. 3) 

E " ' § = T ' v = 2 (B4) 

The expression for L 0 /pZ thus becomes 



Thus from equation (Bl) 

1-E 2 = (1-E)(1 + E) = 2(1 - ■— (B 6 ) 
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For E near unity 



It is clear that E approaches unity as the path becomes vertical. 

At the other limit assume that is nearly equal to V c so that 


V< 


9 


V r 


= 1 + e 


(B8) 


where e is a small quantity. Let us rewrite equation (B7) in the form 


E 2 = 


tp_ 

y c 2 


- 1 + 


m 


(B9) 


or 


E 2 = 


2e 1 + 


12 / . \2 


.v (1+2e+e2) 


(BIO) 


For V C p = V c or very close to it, we thus have 


IV, 


9 


+ 4e 2 
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APPENDIX C 

CHANGE IN PERIGEE ALTITUDE DUE TO OBLATENESS 
GRAVITATIONAL FIELD FOR EQUATORIAL ORBITS 


The specific question considered here corcerns a satellite in a 
circular orbit without oblateness forces at time t = 0: When the oblate- 

ness gravitational field is suddenly "turned cn," what happens to the 
perigee radius, r p ? Let r a be the radius oi the circular orbit without 
oblateness forces, which becomes the apogee radius as soon as oblateness 
forces are turned on. From the previous appendix it is known that the 
total energy is related to the length 2 l of the major axis (without 
oblateness) by equation (b4). 


E = - 


KM 

2Z 


(Cl) 


A change in energy A E is such that 


AE A l 

E " l 


(C2) 


Now for an equatorial orbit the energy balance with . V 9 = V c and with 
no oblateness is such that the orbit is circular. With oblateness the 
energy balance requires that V<p = Vco for a circular orbit. Thus if 
Vcp = V c vith oblateness, there is a kinetic energy deficiency of 
(V co 2 ” V c 2 )/2 per unit mass for a circular orbit. This energy deficiency 
causes a change in length of the major axis L2X of r a - rp since the 
orbit has the same apogee radius with or without oblateness. Thus 


with 


r a “ r p 


2ZAE 

E 



E ~ - 





> 


equation (C3) becomes 



(C3) 


(c4) 


(C5) 


( r a “ r p ) ~ 4P| 
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TABLE I.- PERCENTAGE EFFECT OF OBLATENESS GRAVITATIONAL FIELD ON RANGE 


Cj>A/m = 1 

V(V 0 

0 

-0.01 

-0.05 

-0.10 

miles 

-4406.562 

3582.562 

948.577 

526.382 

A 4j , miles 

20,060.089 

2614.743 

545.045 

335.566 

s 4 , miles 

19,279.667 

4201. 140 

1089.919 

624.226 

t, sec 

4355 

1199 

538 

426 

(^2 " ^l)/s 4 

-.330 

-.042 

-.003 

-.001 

(^2 ™ -^1 )/®4 

-.365 

-.067 

-.001 

.000 

(+4 ‘ + 3 )/ B 4 

-.366 

-.042 

-.003 

-.003 

(A 4 - A 3 )/ s 4 

-.334 

-.070 

-.001 

-.005 

C D A/m = 10 

\[r 4 , miles 

3344.785 

2527.760 

751.924 

427 • 021 

A 4 , miles 

10,539.586 

1689.927 

640.422 

480.512 

s 4 , miles 

8975.380 

2974.090 

984.494 

642. 061 

t , sec 

2801 

1594 

1158 

1068 

(+2 “ t;L)/s 4 

.090 

-.023 

-.002 

-.001 

(A 2 - A 1 )/s 4 

-.095 

-.018 

-.001 

.000 

(t 4 - t 3 )/s 4 

.094 

-.024 

-.002 

-.001 

(A 4 - A 3 )/s 4 

-.091 

-.019 

-.001 

.000 



















TABLE II.- PERCENTAGE EFFECT OF ATMOS PHERCC ROTATION ON RANGE 


CdA/hi = 1 


r 0 /(V 0 

0 

-0.01 

4r 4 , miles 

-4406.562 

3582.562 

A 4 , miles 

20,060.089 

2614.743 

s 4 , miles 

19,279.667 

4201. 140 

t y sec 

4355 

1199 

(i 3 " ti)/s 4 

.o 4 o 

.006 

(A 3 - Ai )/s 4 

.033 

.034 

(i 4 " ^s)/ S 4 

.004 

.005 

(A 4 " A 2 )/ s 4 

.064 

.032 

O 

9 

II 

H 

O 

miles 

3344.785 

2527.760 

A 4 , miles 

10,539.586 

1689.927 

s 4 , miles 

8975 . 380 

2974.090 

t , sec 

2801 

1594 

(+3 " ^l)/ S 4 

-.017 

.007 

(A 3 - A 1 )/s 4 

.045 

.103 

(+4 " ^ 2 )/ S 4 

-.013 

.006 

(A 4 - A 2 )/ s 4 

.049 

.103 


4-0 . 577 

45.045 


1089.919 


751.924 

640.422 

984.494 

1158 

.003 

.291 

.003 

.291 


526.382 

335.566 

624.226 

426 

.oo 4 

.148 

.002 

.143 


427.021 
480.512 
642. 061 
1068 
.002 
.438 
.002 
.438 
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TABLE III.- CALCULATED QUANTITIES FOR EQUATORIAL ORBIT 


t, 

sec 

h, 

ft 

y 

ft/ sec 

-r, 

ft/sec 

A, 

radians 

v c , 

ft/sec 

yv c 

0 

422,400 

25,6780860 

0 

0 

25,678.86 

1.000000 

6 06 

421,753 

25,674.298 

3.899 

.7288 

25,679.25 

.999807 

1669 

407,830 

25,677.980 

24.045 

2.0076 

25,687.62 

. 999624 

1955 

400,009 

25,682.181 

30.820 

2.3519 

25,692.33 

.999605 

2442 

381,881 

25,689.336 

44.438 

2.9387 

25,703.27 

.999458 

2588 

375,008 

25,690.380 

49.84o 

3.1147 

25,707.41 

.999338 

2916 

355,947 

25,685.460 

68.485 

3.5104 

25,718.92 

.998699 

2998 

350,039 

25,680.795 

75.852 

3.6094 

25,722.50 

.998378 

3200 

332,175 

25,651.831 

104 . 165 

3.8532 

25,740.28 

.996834 

3400 

306,036 

25,543.869 

167.173 

4.0943 

25,749.14 

.991989 

3500 

286,112 

25,315.656 

239.939 

4.2142 

25,761.21 

. 982704 

3665 

220,154 

22,022.188 

696.182 

4.4o43 

25,801.36 

.853528 

3696 

195,186 

19,039.298 

923.976 

4.4346 

25,816.60 

.737483 

3721 

169,602 

15,017.882 

1,116.853 

4.4550 

25,832.27 

.581361 

3746 

l4o,388 

9,474.258 

1,183.196 

4.4696 

25,850.16 

.366507 

3755 

129,886 

7,368.195 

1,144.695 

4.4732 

25,856.61 

. 284964 

3764 

119,859 

5,464.839 

1,079-841 

4.4759 

25,862.77 

.211301 

3773 

110,478 

3,896.195 

1,003.839 

4.4779 

25,868.54 

.150615 

3784 

99,935 

2,479.298 

915.170 

4.4796 

24,875.03 

.095818 

3795 

90,295 

1,529.855 

840.117 

4.48o6 

25,880.96 

.059111 

3805 

82,195 

960.579 

781.004 

4.4812 

25,885.94 

.037108 

3821 

70,436 

422.641 

688.208 

4.4817 

25,893.28 

.016315 

3837 

60,268 

159.390 

580.798 

4.4819 

25,899.46 

.006154 

3857 

49,986 

35.608 

452.853 

4.4820 

25,905.81 

. 001374 

3882 

40,103 

3.564 

3^7.173 

4 . 4820 

25,911.91 

.000138 

3915 

30,011 

0O88 

274.510 

4.4820 

25,918.16 

. 000003 

3934 

25,041 

.008 

249.848 

4.4820 

25,921.23 

0 

3955 

20,023 

0 

228.336 

4.4820 

25,924.33 

0 

3978 

14,998 

0 

1 

209.522 

4.4820 

25,927.44 

0 

4003 

9,973 

0 

192.983 

4.4820 

25,930.56 

0 

4030 

4,965 

0 

178.392 

4.4820 

25,933.65 

0 

4058 

151 

0 

165.883 

4.4820 

| 25,936.64 

0 
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Figure 1.- Coordinate systems and notation. 
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(a) Cjyk/m - 1; cases 1 and 3J jjl = 0 • 

Figure 3*" Periodic variations of altitude and eccentricity for 
equatorial orbits; h G = 120 statute miles. 
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(b) Cj)A/m = 10; cases 1 and 3> n = 0. 


Figure 3 • “ Continued 









m = 10; cases 2 and k; |i / 0 


Figure 3 *~ Concluded* 





Eccentricity, ExIO 6 Altitude 


Figure 




2t 

(a) Ci)A/m = 1; cases 1 and 3; \i = 0. 

4.- Secular trends in altitude and eccentricity for equatorial 
orbits; h Q = 120 statute miles. 












m = 10; cases 2 and 4; \i ^ 0 
Figure 4.- Concluded. 
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(a) C D A/m = 10; case 1; n e = 0, q = 0. 

Figure 5*~ Eccentricity variation for equatorial orbits; h D = 120 statute 

miles . 
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Figure 5-" Concluded 
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(a) Position of line of apsides, 6, as polar angle. 

Figure 6.- Position of line of apsides; GpA/m = 10 , h Q = 120 statute 

miles, case l;^ e =0jM = 0. 
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Longitude X or latitude 
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Time, t, seconds 


(b) Latitude and longitude. 


Figure 7." Continued. 



Eccentricity, 


74 



( c ) Eccentric ity « 


Figure 7-" Concluded 
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(c) CdA/iti = 1; r 0 = -0.l(V cp ) o . 


Figure 9 *~ Concluded, 
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Figure 10.- Variation of flight path angle for various observers and conditions of flight; 

a Q = 65°, C-qA/th = 1, h Q = 80 statute miles. 



Altitude, hxlO , feet 


82 


lillfSillS 


^§1 


■■■■■■■■ jililllillil! I 
IBM gjgj g— | ffl ™»|j. 


MMMMPMMMMHBPM 

w ligMMMMMMP 

ImiStliSttml ijiflilj ffiRillH aUjliffl gig HM {ggim UNI 

Imipii mttSHlI sHlIitlii SiRHln BIHfin sIHlIII HtfilHfi nn! 


iissiii mmm 

HH bei ntfigiff 1 jftii twi CTHl 8Hil BHtiMlir 

■ ffiimiH ill! iiHIRIIM, 

■■■■■■I HKiPMiHI 


iHlilllM ■ 
BilffliiiliiiiipPililil'Pi! 

tiixwi imhmhi isiitziui ]uii MiHixin uni txuimii titi; 


mmm liiiiippplpi 




Present method i 

Method of appendix A l 
Large time interval j. 

Small time interval 1 




mmm 


2000 


t, sec 


(a) Altitude. 


Figure 11.- Comparison of approximate analytical trajectory and numerical 
trajectory during initial phase; a 0 = 0°, C 3 A/m =1, 0 e = 0, p = 0, 
h 0 = 80 statute miles. 
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(b) Radial velocity. 


Figure 11.- Concluded. 












